package MathLib.MyMatrices;

public class SLAR extends CMatrix{
    protected static int maxIterNum = 1000000; 
	protected double[] b;
	protected double det;
    protected double eps;
    private int CountOfReplaces;
    protected double[][] U;
    protected double[][] L;
    
    protected double[] a;
    protected double[] c;
    protected double[] f;
    
    
    public SLAR(){
    	super();
    	this.b = null;
    	this.det = 0;
    	this.eps = 0.00000001;
    	this.CountOfReplaces = 0;
    } 
    
    public SLAR(int n,double eps, double[] b){
    	super(n);
    	this.b = (double[]) b.clone();
    	this.det = 0;
    	this.eps = eps;
    	this.CountOfReplaces = 0;
    }
    
    public SLAR(int n,double eps){
    	super(n);
    	this.b = null;
    	this.det = 0;
    	this.eps = eps;
    	this.CountOfReplaces = 0;
    }
    
    public SLAR(CMatrix A,double eps) throws MatrixEception{
    	super(A);
    	if (A.n!=A.m) throw new MatrixEception(3);    	
    	this.b = null;
    	this.det = 0;
    	this.eps = eps;
    	this.CountOfReplaces = 0;
    }
    
    public SLAR(CMatrix A,double eps, double[] b) throws MatrixEception{
    	super(A);
    	if (A.n!=A.m) throw new MatrixEception(3);    	
    	this.b = (double[]) b.clone();
    	this.det = 0;
    	this.eps = eps;
    	this.CountOfReplaces = 0;
    }
    
    public SLAR(CMatrix A,double[] b) throws MatrixEception{
    	super(A);
    	if (A.n!=A.m) throw new MatrixEception(3);    	
    	this.b = (double[]) b.clone();
    	this.det = 0;
    	this.eps = 0.00000001;
    	this.CountOfReplaces = 0;
    }
    
    public SLAR(double[] a,double[] b,double[] c,double[] f){
    	super();
    	this.a = (double[]) a.clone();
    	this.f = (double[]) f.clone();
    	this.c = (double[]) c.clone();
    	this.b = (double[]) b.clone();
    	this.det = 0;
    	this.n = this.m = c.length;
    	this.eps = 0.00000001;
    	this.CountOfReplaces = 0;
    	
    	
    }
    
	
	private  void ChangeOfRows(int i) throws MatrixEception{
		int j;
		for(j = i;(j<n)||(Math.abs(this.Matrix[j][i])>eps);j++){	}
		if(this.Matrix[j][i]!=0){
			double temp;
			for(int k = i;k<n;k++){
				temp = this.Matrix[i][k];
				this.Matrix[i][k] = this.Matrix[j][k];
				this.Matrix[j][k] = temp;
			}		
			temp = b[i];
			b[i]= b[j];
			b[j] = temp;
			CountOfReplaces++;
		}else{
			throw new MatrixEception(1);			
		}
	}
	
	private void mainElement(int i)throws MatrixEception{
		int j,TempPos=i;
		double temp = this.Matrix[i][i];
		for(j=i+1;j<n;j++){
			if(this.Matrix[j][i]>temp){temp = this.Matrix[j][i]; TempPos = j;}
		}
		if (Math.abs(temp)<eps) throw new MatrixEception(1);
		for(j = i;j<n;j++){
			temp = this.Matrix[i][j];
			this.Matrix[i][j] = this.Matrix[TempPos][j];
			this.Matrix[TempPos][j] = temp;
		}		
		temp = b[i];
		b[i]= b[TempPos];
		b[TempPos] = temp;
		CountOfReplaces++;
		
	}
	
	private void mainElement(int i,double[][] B)throws MatrixEception{
		int j,TempPos=i;
		double temp = this.Matrix[i][i];
		for(j=i+1;j<n;j++){
			if(this.Matrix[j][i]>temp){temp = this.Matrix[j][i]; TempPos = j;}
		}
		if (Math.abs(temp)<eps) throw new MatrixEception(1);
		for(j = i;j<n;j++){
			temp = this.Matrix[i][j];
			this.Matrix[i][j] = this.Matrix[TempPos][j];
			this.Matrix[TempPos][j] = temp;
		}		
		
		for(j = 0;j<n;j++){
			temp = B[i][j];
			B[i][j] = B[TempPos][j];
			B[TempPos][j] = temp;
		}
		CountOfReplaces++;
		
	}
    
    public double[] Gauss() throws MatrixEception{
    	if (b == null){throw new MatrixEception(2);}
    	for (int k = 0; k < n-1; k++) {
    		if (Math.abs(this.Matrix[k][k])<eps){
    			ChangeOfRows(k);    			
    		}
    		for(int i=k+1;i<n;i++){
    			
    			double koef = -(this.Matrix[i][k]/this.Matrix[k][k]);
    			for(int j = k+1;j<n;j++){
    				this.Matrix[i][j] += koef*this.Matrix[k][j];
    			}
    			b[i] += koef*b[k];
    		}
		}
        double[] res = new double[n];
        int k = 0;
        double temp;
    	for (int i = n-1; i >= 0; i--) {
    		temp = 0;
    		for(k=i+1;k<n;k++){
    			temp += res[k]*this.Matrix[i][k];
    		}
    		res[i] = (b[i] - temp)/this.Matrix[i][i];
		}
    	
        return res;
    }
    
    public double[] ModefiedGauss() throws MatrixEception{
    	if (b == null){throw new MatrixEception(2);}
    	for (int k = 0; k < n-1; k++) {
    		mainElement(k);
    		for(int i=k+1;i<n;i++){
    			
    			double koef = -(this.Matrix[i][k]/this.Matrix[k][k]);
    			for(int j = k+1;j<n;j++){
    				this.Matrix[i][j] += koef*this.Matrix[k][j];
    			}
    			b[i] += koef*b[k];
    		}
		}
        double[] res = new double[n];
        int k = 0;
        double temp;
    	for (int i = n-1; i >= 0; i--) {
    		temp = 0;
    		for(k=i+1;k<n;k++){
    			temp += res[k]*this.Matrix[i][k];
    		}
    		res[i] = (b[i] - temp)/this.Matrix[i][i];
		}
    	
        return res;
    }
    
    public CMatrix Tranverse_Gauss() throws MatrixEception{
    	//if (b == null){throw new MatrixEception(2);}
    	double[][] B = new double[n][n]; //// ��������� �� ������ �������� ���� ���������;
    	for(int i = 0;i<n;i++){
    		B[i][i] = 1;    		
    	}
    	
    	for (int k = 0; k < n-1; k++) {
    		mainElement(k,B);
    		for(int i=k+1;i<n;i++){
    			
    			double koef = -(this.Matrix[i][k]/this.Matrix[k][k]);
    			for(int j = k+1;j<n;j++){
    				this.Matrix[i][j] += koef*this.Matrix[k][j];
    			}
    			for(int j = 0;j<n;j++){
    				B[i][j] += koef*B[k][j];
    			}
    		}
		}
    	CMatrix A = new CMatrix(n);
    	double[] res = new double[n];
    	for(int j = 0;j<n;j++){
    		int k = 0;
        	double temp;
    		for (int i = n-1; i >= 0; i--) {
    			temp = 0;
    			for(k=i+1;k<n;k++){
    				temp += res[k]*this.Matrix[i][k];
    			}
    			res[i] = (B[i][j] - temp)/this.Matrix[i][i];
			}
    		//A.Matrix[j] = res;
    		for(int l=0;l<n;l++){
    			A.Matrix[l][j] = res[l];    			
    		}
    	}
    	return A;
    }
    
    public CMatrix Tranverse_LU() throws MatrixEception{
    	if((L==null)||(U==null)){LUFactorization();}
    	double[] y = new double[n];
    	double[] x = new double[n];
    	CMatrix A = new CMatrix(n);
    	double temp;
    	for(int q = 0;q<n;q++){  	    	
    			for(int i=0;i<q;i++){
    				temp=0;
    				for (int k = 0; k <= i-1; k++) {
    	    			temp+=L[i][k]*y[k];
    	    		}
    				y[i] = (-temp);
    			}
    			
    				temp=0;
    				for (int k = 0; k <= q-1; k++) {
    					temp+=L[q][k]*y[k];
    				}
    				y[q] = (1-temp);
    			
    			
    			for(int i=q+1;i<n;i++){
    				temp=0;
    				for (int k = 0; k <= i-1; k++) {
    	    			temp+=L[i][k]*y[k];
    	    		}
    				y[i] = (-temp);
    			}
    			
    			for(int i=n-1;i>=0;i--){
    				temp=0;
    				for (int k = i+1; k < n; k++) {
    					temp+=U[i][k-i]*x[k];
    				}
    				x[i] = (y[i]-temp)/U[i][0];
    			}  	
    		for(int l=0;l<n;l++){
        		A.Matrix[l][q] = x[l];    			
        	}
    	}
    	
    	
    	return A;
    }
    
    public double[] SqareRootMethod() throws MatrixEception{
    	if(!IsSumetr()){new MatrixEception(3);}
    	U = new double[n][n];
    	
    	for (int i = 0; i < n; i++) {
    		double temp = 0;
    		for (int k = 0; k < i; k++) {
    			temp+=U[k][i]*U[k][i];
    		}
    		if(Math.abs(Matrix[i][i]-temp)<eps) throw new MatrixEception(4);
    			U[i][i] = Math.sqrt(Matrix[i][i]-temp);
    		
    		for (int j = 1; j < n; j++) {
    			temp = 0;
    			for (int k = 0; k < i; k++) {
        			temp+=U[k][i]*U[k][j];
        		}
    			try{
    			U[i][j] = (Matrix[i][j]-temp)/U[i][i];
    			}catch(Exception e){
    				throw new MatrixEception(5);   				
    			}
    		}
		}    	
    	
    	double[] y = new double[n];
    		for(int i=0;i<n;i++){
    			double temp=0;
    			for (int k = 0; k < i-1; k++) {
        			temp+=U[k][i]*y[k];
        		}
    			y[i] = (b[i]-temp)/U[i][i];
    		}
    		
    	double[] x = new double[n];
    		for(int i=n-1;i>=0;i--){
    			double temp=0;
    			for (int k = i+1; k < n; k++) {
    				temp+=U[i][k]*x[k];
    			}
    			x[i] = (y[i]-temp)/U[i][i];
    		}
    		
    	
    	
    	return x;
    }
    
    public void LUFactorization() throws MatrixEception{
    	U = new double[n][];
    	L = new double[n][];
    		for(int i=1;i<=n;i++){
    			U[i-1] = new double[n-i+1];
    			L[i-1] = new double[i];
    			U[0][i-1] = Matrix[0][i-1];
    			L[i-1][i-1] = 1;
    		}
    		for(int i=0;i<n;i++){

//    			L[][]
    			
    			for(int j=0;j<i;j++){
    				double temp = 0;
    				for (int k = 0;k<i-1;k++){
    					temp+=L[i][k]*U[k][j-k];
    				}
    				if(Math.abs(U[j][0])<eps){throw new MatrixEception(6);}
    				L[i][j] = (Matrix[i][j] - temp)/U[j][0];
    			}
//    			u[][]
    			for(int j=i;j<n;j++){
    				double temp = 0;
    				for (int k = 0;k<=i-1;k++){
    					temp+=L[i][k]*U[k][j-k];
    				}
    				U[i][j-i] = Matrix[i][j] - temp;
    			}
    		}
    }
    
    public double[] SolveUsingLU(){
    	double[] y = new double[n];
		for(int i=0;i<n;i++){
			double temp=0;
			for (int k = 0; k <= i-1; k++) {
    			temp+=L[i][k]*y[k];
    		}
			y[i] = (b[i]-temp);
		}
		
		double[] x = new double[n];
		for(int i=n-1;i>=0;i--){
			double temp=0;
			for (int k = i+1; k < n; k++) {
				temp+=U[i][k-i]*x[k];
			}
			x[i] = (y[i]-temp)/U[i][0];
		}
    	return x;
    }
    
    public double[] LeftMarching() {
    	
    	double[] Alfa = new double[n];
    	double[] Beta = new double[n];
    	
    	Alfa[n-1]=(-c[n-1])/a[n-1];
    	Beta[n-1]=(f[n-1])/a[n-1];

    	
    	for (int i=(n-2); i>(-1); i--)
    	   {
    	       Alfa[i]=(-b[i+1]-Alfa[i+1]*c[i+1])/(Alfa[i+1]*a[i+1]);
    	       Beta[i]=(f[i+1]-Beta[i+1]*c[i+1]-Beta[i+1]*a[i+1]*Alfa[i])/(a[i+1]);
    	   }
    	double[] x = new double[n];

    	   x[0]=(b[0]*Beta[0]+f[0]*Alfa[0])/(c[0]*Alfa[0]+b[0]);

    	   for (int i=2; i<(n+1); i++)
    	   {
    	       x[i-1]=(x[i-2]-Beta[i-2])/Alfa[i-2];
    	   }
    	
		return x;
	}
    
    public double[] RightMarching() {
    	//int parameter = 1;
    	double[] Alfa = new double[n];
    	double[] Beta = new double[n];
    	
    	Alfa[0]=(-b[0])/c[0];
    	Beta[0]=f[0]/c[0];
    	   
    	for (int i=1; i<n; i++)
    	   {
    	       Alfa[i]=(-b[i])/(c[i]+Alfa[i-1]*a[i]);
    	       Beta[i]=(f[i]-Beta[i-1]*a[i])/(c[i]+Alfa[i-1]*a[i]);
    	   }
    	double[] x = new double[n];
    	   x[n-1]=(f[n-1]-Beta[n-2]*a[n-1])/(c[n-1]+Alfa[n-2]*a[n-1]);
    	   
    	   for (int i=(n-2); i>=(0); i--)
    	   {
    	       x[i]=Alfa[i]*x[i+1]+Beta[i];
    	   }
    	
		return x;
	}
    
    private double[] LeftMarching(int param,double[] x) {
    	
    	double[] Alfa = new double[n];
    	double[] Beta = new double[n];
    	
    	Alfa[n-1]=(-c[n-1])/a[n-1];
    	Beta[n-1]=(f[n-1])/a[n-1];

    	
    	for (int i=(n-2); i>(-1); i--)
    	   {
    	       Alfa[i]=(-b[i+1]-Alfa[i+1]*c[i+1])/(Alfa[i+1]*a[i+1]);
    	       Beta[i]=(f[i+1]-Beta[i+1]*c[i+1]-Beta[i+1]*a[i+1]*Alfa[i])/(a[i+1]);
    	   }

    	   x[0]=(b[0]*Beta[0]+f[0]*Alfa[0])/(c[0]*Alfa[0]+b[0]);

    	   for (int i=2; i<(param); i++)
    	   {
    	       x[i-1]=(x[i-2]-Beta[i-2])/Alfa[i-2];
    	   }
    	
		return x;
	}
    
    private double[] RightMarching(int param,double[] x) {
    	//int parameter = 1;
    	double[] Alfa = new double[n];
    	double[] Beta = new double[n];
    	
    	Alfa[0]=(-b[0])/c[0];
    	Beta[0]=f[0]/c[0];
    	   
    	for (int i=1; i<n; i++)
    	   {
    	       Alfa[i]=(-b[i])/(c[i]+Alfa[i-1]*a[i]);
    	       Beta[i]=(f[i]-Beta[i-1]*a[i])/(c[i]+Alfa[i-1]*a[i]);
    	   }
    	   x[n-1]=(f[n-1]-Beta[n-2]*a[n-1])/(c[n-1]+Alfa[n-2]*a[n-1]);
    	   
    	   for (int i=(n-2); i>=(param-1); i--)
    	   {
    	       x[i]=Alfa[i]*x[i+1]+Beta[i];
    	   }
    	
		return x;
	}
    
    public double[] MeetMarching(int param){
    	
    	double[] x = new double[n];
    	LeftMarching(param, x);
    	RightMarching(param, x);
    	return x;
    }
    
    public double[] RotationMethod() throws MatrixEception{
    	if (b == null){throw new MatrixEception(2);}
    	for (int k = 0; k < n-1; k++) {
    		if (Math.abs(this.Matrix[k][k])<eps){
    			ChangeOfRows(k);    			
    		}
    		
    		for(int i=k+1;i<n;i++){
    			if(Math.abs(this.Matrix[i][k])<eps){continue;}
    			double c = Matrix[k][k]/(Math.sqrt(Matrix[k][k]*Matrix[k][k]+Matrix[i][k]*Matrix[i][k]));
    			double s = Matrix[i][k]/(Math.sqrt(Matrix[k][k]*Matrix[k][k]+Matrix[i][k]*Matrix[i][k]));
    			double temp,tmp;
    			for(int j = k;j<n;j++){
    				
    				temp = Matrix[k][j];
    				tmp = Matrix[i][j];
    				Matrix[k][j] = c*temp+s*tmp;
    				Matrix[i][j] = -s*temp+c*tmp;
    				
    			}
    			temp = b[k];
				tmp = b[i];
				b[k] = c*temp+s*tmp;
				b[i] = -s*temp+c*tmp;
    			
    		}
		}
    	double[] res = new double[n];
        int k = 0;
        double temp;
    	for (int i = n-1; i >= 0; i--) {
    		temp = 0;
    		for(k=i+1;k<n;k++){
    			temp += res[k]*this.Matrix[i][k];
    		}
    		res[i] = (b[i] - temp)/this.Matrix[i][i];
		}
    	
        return res;
    }
    
    public double[] JakobiMethod(double[] x0,double epss){
    	int iter = 0;
    	boolean key = true;
    	double[] res = JakobiM( x0, epss);
    	while(key&&iter!=maxIterNum){
    		iter++;
    		key = false;
    		for (int i=0; i<n; i++){
    			double tmp=0.0;
    			for (int j=0; j<n; j++){
    				tmp+=this.Matrix[i][j]*res[j];
    			}
    			
    			if(Math.abs(tmp-b[i])>epss){
    				key=true;
    				res = JakobiM( res, epss);
    				break;
    			}
    		}
    	}
    	return res;
    }
    
    public double[] JakobiM(double[] x0,double epss) {
    	double[] d = new double[n];
		CMatrix C = new CMatrix(Matrix);
		for(int i=0;i<n;i++){
			C.set(i, i, 0);
			d[i] = b[i]/Matrix[i][i];
		}
		
		for(int i =0;i<n;i++){
			for(int j=0;j<n;j++){
				C.set(i, j, -C.get(i, j)/Matrix[i][i]);
			}
		}
		
		boolean key = true;
		double[] x = new double[n];
		while(key){
			for(int i =0;i<n;i++){
				double tmp = 0.0;
				for(int j=0;j<n;j++){
					tmp+=C.get(i, j)*x0[j];
				}
				x[i] = tmp+d[i];
			}
			key = false;
			for(int j=0;j<n;j++){
				if(Math.abs(x[j]-x0[j])>epss){
					key = true;
					x0 = x;
				}
			}
		}
		
		return x;
	}
    
    public double[] ZejdelMethod(double[] x0,double epss) throws MatrixEception{
    	int iter = 0;
    	boolean key = true;
    	double[] res = ZejdelM( x0, epss);
    	while(key&&iter!=maxIterNum){
    		iter++;
    		key = false;
    		for (int i=0; i<n; i++){
    			double tmp=0.0;
    			for (int j=0; j<n; j++){
    				tmp+=this.Matrix[i][j]*res[j];
    			}
    			
    			if(Math.abs(tmp-b[i])>epss){
    				key=true;
    				res = ZejdelM( res, epss);
    				break;
    			}
    		}
    	}
    	return res;
    }

    public double[] ZejdelM(double[] x0,double epss) throws MatrixEception {
    	double[] d = new double[n];
    	CMatrix C = new CMatrix(Matrix);
    	CMatrix R = new CMatrix(Matrix);

		for(int i =0;i<n;i++){
			
			for(int j=0;j<=i;j++){
				R.set(i, j, 0);
			}
			for(int j=i+1;j<n;j++){
				C.set(i, j, 0);
			}
		}
		
		SLAR temp = new SLAR(n,0.000001);
		temp.Matrix = C.Matrix;
		C = temp.Tranverse_LU();
		//C.Matrix = temp.Matrix;
		
		for(int i =0;i<n;i++){
			double z=0.0;
			for(int j=0;j<n;j++){
				z+=C.Matrix[i][j]*b[j];
			}
			d[i] = z;
		}
		
		C.Multiply(R);
		
		for(int i =0;i<n;i++){
			

			for(int j=0;j<i+1;j++){
				C.Matrix[i][j]*=-1;
			}
		}
		
		boolean key = true;
		double[] x = new double[n];
		while(key){
			for(int i =0;i<n;i++){
				double tmp = 0.0;
				for(int j=0;j<i;j++){
					tmp+=C.get(i, j)*x[j];
				}
				for(int j=i;j<n;j++){
					tmp+=C.get(i, j)*x0[j];
				}
				x[i] = tmp+d[i];
			}
			key = false;
			for(int j=0;j<n;j++){
				if(Math.abs(x[j]-x0[j])>epss){
					key = true;
					x0 = x;
				}
			}
		}
		
		return x;
	}
    
    public double[] RelacsationMethod(double[] x0,double epss,double w) throws MatrixEception{
    	int iter = 0;
    	boolean key = true;
    	double[] res = RelacsationM( x0, epss,w);
    	while(key&&iter!=maxIterNum){
    		iter++;
    		key = false;
    		for (int i=0; i<n; i++){
    			double tmp=0.0;
    			for (int j=0; j<n; j++){
    				tmp+=this.Matrix[i][j]*res[j];
    			}
    			
    			if(Math.abs(tmp-b[i])>epss){
    				key=true;
    				res = RelacsationM( res, epss,w);
    				break;
    			}
    		}
    	}
    	return res;
    }
    
    public double[] RelacsationM(double[] x0,double epss,double w) throws MatrixEception {
    	double[] d = new double[n];
    	CMatrix C = new CMatrix(Matrix);
    	CMatrix R = new CMatrix(Matrix);

		for(int i =0;i<n;i++){
			
			for(int j=0;j<=i;j++){
				R.set(i, j, 0);
			}
			for(int j=i+1;j<n;j++){
				C.set(i, j, 0);
			}
		}
		
		SLAR temp = new SLAR(n,0.000001);
		temp.Matrix = C.Matrix;
		C = temp.Tranverse_LU();
		
		for(int i =0;i<n;i++){
			double z=0.0;
			for(int j=0;j<n;j++){
				z+=C.Matrix[i][j]*b[j];
			}
			d[i] = z;
		}
		
		C.Multiply(R);
		
		for(int i =0;i<n;i++){

			for(int j=0;j<i+1;j++){
				C.Matrix[i][j]*=-1;
			}
		}
		
		boolean key = true;
		double[] x = new double[n];
		double[] z = new double[n];
		while(key){
			for(int i =0;i<n;i++){
				double tmp = 0.0;
				for(int j=0;j<i;j++){
					tmp+=C.get(i, j)*z[j];
				}
				for(int j=i;j<n;j++){
					tmp+=C.get(i, j)*x0[j];
				}
				z[i] = tmp+d[i];
			}
			for(int i =0;i<n;i++){
				z[i]-=x0[i];
				z[i]*=w;
			}
			for(int i =0;i<n;i++){
				x[i]=x0[i]+z[i];
			}
			key = false;
			for(int j=0;j<n;j++){
				if(Math.abs(x[j]-x0[j])>epss){
					key = true;
					x0 = x;
				}
			}
		}
		
		return x;
	}
    
    public double[] MSG(double[] x0,double epss) throws MatrixEception{
    	double[] r = new double[n];
    	double[] p = new double[n];
    	double[] x = new double[n];
    	double[] Ap = new double[n];
    	
    	for(int i =0;i<n;i++){
			double tmp = 0.0;
			for(int j=0;j<n;j++){
				tmp+=Matrix[i][j]*x0[j];
			}
			r[i] = b[i]-tmp;
		}
    	p=r;
    	boolean key = true;
    	while(key){
    		double alpha=0.0;
    		
    		for(int i =0;i<n;i++){
    			alpha+=r[i]*p[i];
    		}
    		for(int i =0;i<n;i++){
    			double tmp = 0.0;
    			for(int j=0;j<n;j++){
    				tmp+=Matrix[i][j]*p[j];
    			}
    			Ap[i] = tmp;
    		}
    		double nuz = 0.0;
    		
    		for(int i =0;i<n;i++){
    			nuz+=p[i]*Ap[i];
    		}
    		alpha/=nuz;
    		
    		for(int i =0;i<n;i++){
    			x[i] = x0[i]+alpha*p[i];
    			r[i] -= alpha*Ap[i];
    		}
    		nuz=0.0;
    		for(int i =0;i<n;i++){
    			nuz+=r[i]*r[i];
    		}
    		if(Math.abs(nuz)<epss){
    			key = false;
    		}else{
    			nuz = 0.0;
    			alpha = 0.0;
    			for(int i =0;i<n;i++){
        			nuz+=p[i]*Ap[i];
        			alpha+=r[i]*Ap[i];
        		}
    			double beta = alpha/nuz;
    			for(int i =0;i<n;i++){
        			p[i] = r[i] - beta*p[i];
        		}
    			x0 = x;
    		}
    	}
    	
    
    	
    	return x;
    }
    //public 
}

/*
 * */




